First we import our tidied salamander morphology & trait dataset.
# Import tidied salamander eye size and trait data
salamanders <- data.frame(read.csv("../Data/Tidy/salamanders_tidy.csv", header=TRUE, na.strings=c("", "NA", " ","<1")))
# Quick look at data structure
str(salamanders)
Next, we import the amphibian tree from Jetz and Pyron (2019) that has been pruned to match the salamander species in this dataset.
# Import pruned phylogeny
caudatatree <- read.nexus("../Data/Tidy/caudata-tree.nex")
# Plot tree
plot.phylo(caudatatree, show.tip.label=TRUE, cex = 0.7)
First, we can take a preliminary look at all of our data so that we can see how it looks and whether we notice any obvious outliers that we might want to go back and check out. Here we plot out our measurements from the right side vs. the left side of each specimen so that we can make sure that we have consistent measurements throughout the dataset. We use the package plotly to create interactive plots - hover over a point to see the name of the species measured.
We expect that measurements of left and right eyes and corneas should be similar within specimens. We can plot out our data to test for this and make sure we don’t see any glaring issues. We can also fit a regression, which we would expect to have a slope near 1 and an intercept near 0 if left and right eye measurements are similar.
#### right eye diameter vs. left eye diameter ####
# Linear model
RLeye.lm <- lm(ED_right_mm ~ ED_left_mm, data = salamanders)
summary(RLeye.lm)
##
## Call:
## lm(formula = ED_right_mm ~ ED_left_mm, data = salamanders)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.56832 -0.13433 -0.01103 0.13191 0.67936
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.016679 0.035501 0.47 0.639
## ED_left_mm 1.009053 0.009215 109.51 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2063 on 391 degrees of freedom
## (45 observations deleted due to missingness)
## Multiple R-squared: 0.9684, Adjusted R-squared: 0.9683
## F-statistic: 1.199e+04 on 1 and 391 DF, p-value: < 2.2e-16
# Plot
RLeye.plot <- ggplot(salamanders, aes(x = ED_left_mm, y = ED_right_mm, text = Genus_Species)) +
geom_point(alpha = 0.9) +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
geom_abline(slope = coef(RLeye.lm)[[2]], intercept = coef(RLeye.lm)[[1]])
# Interactive plot
ggplotly(RLeye.plot)
As we would expect, measurements of the left and right eye are highly correlated, with a slope close to 1 and an intercept close to 0. We can do the same test for the corneal measurements.
#### right cornea diameter vs. left cornea diameter ####
# Linear model
RLcornea.lm <- lm(CD_right_mm ~ CD_left_mm, data = salamanders)
summary(RLcornea.lm)
##
## Call:
## lm(formula = CD_right_mm ~ CD_left_mm, data = salamanders)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.70444 -0.11965 -0.01275 0.14285 0.79425
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.11282 0.03227 3.497 0.00052 ***
## CD_left_mm 0.97073 0.01163 83.458 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2149 on 434 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.9413, Adjusted R-squared: 0.9412
## F-statistic: 6965 on 1 and 434 DF, p-value: < 2.2e-16
# Plot
RLcor.plot <- ggplot(salamanders, aes(x = CD_left_mm, y = CD_right_mm, text = Genus_Species)) +
geom_point(alpha = 0.9) +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
geom_abline(slope = coef(RLcornea.lm)[[2]], intercept = coef(RLcornea.lm)[[1]])
# Interactive plot
ggplotly(RLcor.plot)
Again, we find that left and right corneal measurements are highly correlated, with a slope close to 1 and an intercept close to 0. Together, these results indicate that measurements were consistent between the left and right eyes of a specimen, which is good! For looking at interspecific allometry, we will calculate species means from these data so that we can use phylogenetic comparative methods.
Here, we will use phylogenetic generalised least-squares (PGLS) regressions to see how eyes are scaling with body size across different species of salamanders with different ecologies. PGLS analyses are similar to normal linear regressions, but also account for the non-independence of data collected from close relatives due to shared evolutionary history.
In these analyses, we will log eye size and body size variables, as this is standard for allometric relationships. First, we will calculate the mean measurements across specimens for each species. We need to make sure in doing this that we omit rows with missing data for a variable of interest. In looking at the data, we can see that we have many more measurements for cornea diameter than for eye diameter.
We then find species means for various morphological parameters. There are some specimens that lack data for eye size but have data for svl, mass, and corneal diameter (these measurements are present for all specimens). However, when eye size measurements are missing, they are consistently missing for all specimens measured within a species. This means we don’t have to worry about missing data for one specimen corrupting the mean across several specimens for eye size (it is either present for all specimens in a species or absent for all specimens in a species). So we can make just one dataframe of species averages, and some species will just lack eye size data.
# Species means for eye size, cornea size, SVL, and rootmass data
av.sal <- salamanders %>%
mutate_if(is.character, as.factor) %>%
group_by(Genus_Species) %>%
summarise(eye_av = mean(eyemean),
cor_av = mean(cormean),
svl_av = mean(SVL_mm),
rootmass_av = mean(rootmass),
n = n())
# Merge data means with other column info, keeping only trait columns
av.sal <- merge(av.sal, salamanders[match(unique(salamanders$Genus_Species), salamanders$Genus_Species), ], by="Genus_Species", all.x = TRUE, all.y = FALSE) %>%
select(Genus_Species, eye_av, cor_av, svl_av, rootmass_av, n, Order, Suborder, Family, Genus, Species, Gill_Presence, Activity_Period, Adult_Habitat, Development, Life_History, Larval_Habitat)
# Remove Eurycea rathbuni and Proteus anguinus (blind salamanders) from dataframe
av.sal <- av.sal[-c(77, 133),]
# Export tidied dataset for analyses
write.csv(av.sal, file = "../Data/Tidy/salamanders_averages.csv")
Here we run a custom color palette for the adult habitat categories to be used for all following figures.
First we define what we want that color palette to be:
# See what the levels of adult habitat are named
levels(as.factor(av.sal$Adult_Habitat))
## [1] "aquatic" "scansorial" "semiaquatic" "subfossorial" "terrestrial"
# Define a colorblind-friendly vector of colors for adult habitat
col_hab <- c("aquatic" = "#0072B2",
"scansorial" = "#009E73",
"semiaquatic" = "#56B4E9",
"subfossorial" = "#CC79A7",
"terrestrial" = "#E69F00")
# Now you can see that each habitat is assigned a hex color when you look at the vector
col_hab
## aquatic scansorial semiaquatic subfossorial terrestrial
## "#0072B2" "#009E73" "#56B4E9" "#CC79A7" "#E69F00"
# Make a vector of point shapes for adult habitat
shape_hab <- c("aquatic" = 25,
"scansorial" = 24,
"semiaquatic" = 15,
"subfossorial" = 16,
"terrestrial" = 18)
#see that each state has a shape now
shape_hab
## aquatic scansorial semiaquatic subfossorial terrestrial
## 25 24 15 16 18
In order to fit a PGLS regression in caper, we first need to make a comparative data object that includes our dataset and our phylogeny. Note that the gilled Pleurodeles waltl will be dropped from this, as there can only be one row for each species.
# Make row names of dataset the species names (so it will match phylogeny tips)
rownames(av.sal) <- av.sal$Genus_Species
# Check that names match in dataframe and tree
name.check(phy = caudatatree, data = av.sal)
# Use caper function to combine phylogeny and data into one object (this function also matches species names in tree and dataset)
sal.comp <- comparative.data(phy = caudatatree, data = av.sal,
names.col = Genus_Species, vcv = TRUE,
na.omit = FALSE, warn.dropped = TRUE)
# Check for dropped tips or dropped species
sal.comp$dropped$tips #phylogeny
sal.comp$dropped$unmatched.rows #dataset
Here we fit PGLS regressions and plot the resulting outputs for the logged eye size and body size variables.
# PGLS eye diameter vs. cube root of mass
pgls_edmass <- pgls(log10(eye_av) ~ log10(rootmass_av),
data = sal.comp,
lambda = "ML", #uses Maximum Likelihood estimate of lambda
param.CI = 0.95)
# Diagnostic plots
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_edmass)
par(mfrow = c(1, 1))
# Likelihood plot for Pagel's lambda. Solid red line indicates estimate for lambda and broken red lines indicate the 95% confidence interval
lambda.edmass <- pgls.profile(pgls_edmass, "lambda")
plot(lambda.edmass)
# Print model output
summary(pgls_edmass)
##
## Call:
## pgls(formula = log10(eye_av) ~ log10(rootmass_av), data = sal.comp,
## lambda = "ML", param.CI = 0.95)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0283788 -0.0042682 0.0000067 0.0040428 0.0180521
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.845
## lower bound : 0.000, p = 6.8369e-11
## upper bound : 1.000, p = 9.5763e-07
## 95.0% CI : (0.654, 0.943)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.392719 0.042518 9.2365 4.441e-16 ***
## log10(rootmass_av) 0.798584 0.043509 18.3546 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.00678 on 137 degrees of freedom
## (16 observations deleted due to missingness)
## Multiple R-squared: 0.7109, Adjusted R-squared: 0.7088
## F-statistic: 336.9 on 1 and 137 DF, p-value: < 2.2e-16
# Quick plot of pgls
plot(log10(eye_av) ~ log10(rootmass_av), data = av.sal)
abline(pgls_edmass)
# Plot average eye diameter (mm) against cube root of mass (g)
plot_pgls_edmass <- ggplot(av.sal, aes(y=eye_av, x=rootmass_av, text=Genus_Species)) +
geom_point(alpha = 0.9, aes(color = Adult_Habitat), size = 5, na.rm = TRUE) +
scale_color_manual(values = col_hab, name = "Adult habitat") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_text(size = 26),
axis.text.y = element_text(size = 26),
axis.title.x = element_text(size = 30, color = "black", margin = margin(t = 10)),
axis.title.y = element_text(size = 30, color = "black", margin = margin(r = 18)),
axis.ticks.length = unit(0.3, "cm"),
panel.border = element_rect(size = 2)
) +
scale_y_log10(name = "Eye diameter (mm)", breaks = c(2, 3, 7)) +
scale_x_log10(name = "Cube root of mass (g)", breaks = c(1, 3, 5), limits = c(0.5, 5.5)) +
geom_abline(slope = coef(pgls_edmass)[[2]], intercept = coef(pgls_edmass)[[1]], linetype = "solid", size = 1.5)
plot(plot_pgls_edmass)
# Save the plot as a rectangle with specified width and height
ggsave("../Plots/plot_pgls_edmass.png", plot = plot_pgls_edmass, width = 9, height = 7.5, units = "in")
# PGLS eye diameter vs. SVL
pgls_edsvl <- pgls(log10(eye_av) ~ log10(svl_av),
data = sal.comp,
lambda = "ML", #uses Maximum Likelihood estimate of lambda
param.CI = 0.95)
# Diagnostic plots
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_edsvl)
par(mfrow = c(1, 1))
# Likelihood plot for Pagel's lambda. Solid red line indicates estimate for lambda and broken red lines indicate the 95% confidence interval
lambda.edsvl <- pgls.profile(pgls_edsvl, "lambda")
plot(lambda.edsvl)
# Print model output
summary(pgls_edsvl)
##
## Call:
## pgls(formula = log10(eye_av) ~ log10(svl_av), data = sal.comp,
## lambda = "ML", param.CI = 0.95)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0206434 -0.0040166 0.0003049 0.0037003 0.0175331
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.779
## lower bound : 0.000, p = 5.7065e-14
## upper bound : 1.000, p = 2.6052e-08
## 95.0% CI : (0.542, 0.918)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.887868 0.082795 -10.724 < 2.2e-16 ***
## log10(svl_av) 0.809726 0.041312 19.600 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.005968 on 137 degrees of freedom
## (16 observations deleted due to missingness)
## Multiple R-squared: 0.7371, Adjusted R-squared: 0.7352
## F-statistic: 384.2 on 1 and 137 DF, p-value: < 2.2e-16
# Quick plot of pgls
plot(log10(eye_av) ~ log10(svl_av), data = av.sal)
abline(pgls_edsvl)
# Plot average eye diameter (mm) against snout-vent length (mm)
plot_pgls_edsvl <- ggplot(av.sal, aes(y=eye_av, x=svl_av, text=Genus_Species)) +
geom_point(alpha = 0.9, aes(color = Adult_Habitat), size = 5, na.rm = TRUE) +
scale_color_manual(values = col_hab, name = "Adult habitat") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_text(size = 26),
axis.text.y = element_text(size = 26),
axis.title.x = element_text(size = 30, color = "black", margin = margin(t = 10)),
axis.title.y = element_text(size = 30, color = "black", margin = margin(r = 18)),
axis.ticks.length = unit(0.3, "cm"),
panel.border = element_rect(size = 2)
) +
scale_y_log10(name = "Eye diameter (mm)", breaks = c(2, 3, 7)) +
scale_x_log10(name = "Snout-vent length (mm)", breaks = c(30, 70, 150), limits = c(25,160)) +
geom_abline(slope = coef(pgls_edsvl)[[2]], intercept = coef(pgls_edsvl)[[1]], linetype = "solid", size = 1.5)
plot(plot_pgls_edsvl)
# Save the plot as a rectangle with specified width and height
ggsave("../Plots/plot_pgls_edsvl.png", plot = plot_pgls_edsvl, width = 9, height = 7.5, units = "in")
# PGLS eye diameter vs. cornea diameter
pgls_edcd <- pgls(log10(cor_av) ~ log10(eye_av),
data = sal.comp,
lambda = "ML", #uses Maximum Likelihood estimate of lambda
param.CI = 0.95)
# Diagnostic plots
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_edcd)
par(mfrow = c(1, 1))
# Likelihood plot for Pagel's lambda. Solid red line indicates estimate for lambda and broken red lines indicate the 95% confidence interval
lambda.edcd <- pgls.profile(pgls_edcd, "lambda")
plot(lambda.edcd)
# Print model output
summary(pgls_edcd)
##
## Call:
## pgls(formula = log10(cor_av) ~ log10(eye_av), data = sal.comp,
## lambda = "ML", param.CI = 0.95)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0078114 -0.0015818 -0.0000318 0.0016618 0.0076100
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.514
## lower bound : 0.000, p = 0.00035263
## upper bound : 1.000, p = 8.8929e-14
## 95.0% CI : (0.132, 0.817)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.141399 0.017157 -8.2416 1.223e-13 ***
## log10(eye_av) 0.980989 0.021445 45.7444 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.002427 on 137 degrees of freedom
## (16 observations deleted due to missingness)
## Multiple R-squared: 0.9386, Adjusted R-squared: 0.9381
## F-statistic: 2093 on 1 and 137 DF, p-value: < 2.2e-16
# Quick plot of pgls
plot(log10(cor_av) ~ log10(eye_av), data = av.sal)
abline(pgls_edcd)
# Plot average eye diameter (mm) against cornea diameter (mm)
plot_pgls_edcd <- ggplot(av.sal, aes(y=cor_av, x=eye_av, text=Genus_Species)) +
geom_point(alpha = 0.9, aes(color = Adult_Habitat), size = 5, na.rm = TRUE) +
scale_color_manual(values = col_hab, name = "Adult habitat") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_text(size = 26),
axis.text.y = element_text(size = 26),
axis.title.x = element_text(size = 30, color = "black", margin = margin(t = 10)),
axis.title.y = element_text(size = 30, color = "black", margin = margin(r = 18)),
axis.ticks.length = unit(0.3, "cm"),
panel.border = element_rect(size = 2)
) +
scale_y_log10(name = "Cornea diameter (mm)") +
scale_x_log10(name = "Eye diameter (mm)", breaks = c(2, 3, 7)) +
geom_abline(slope = coef(pgls_edcd)[[2]], intercept = coef(pgls_edcd)[[1]], linetype = "solid", size = 1.5)
plot(plot_pgls_edcd)
# Save the plot as a rectangle with specified width and height
ggsave("../Plots/plot_pgls_edcd.png", plot = plot_pgls_edcd, width = 9, height = 7.5, units = "in")
# PGLS cornea diameter vs. cube root of mass
pgls_cdmass <- pgls(log10(cor_av) ~ log10(rootmass_av),
data = sal.comp,
lambda = "ML", #uses Maximum Likelihood estimate of lambda
param.CI = 0.95)
# Diagnostic plots
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_cdmass)
par(mfrow = c(1, 1))
# Likelihood plot for Pagel's lambda. Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval
lambda.cdmass <- pgls.profile(pgls_cdmass, "lambda")
plot(lambda.cdmass)
# Print model output
summary(pgls_cdmass)
##
## Call:
## pgls(formula = log10(cor_av) ~ log10(rootmass_av), data = sal.comp,
## lambda = "ML", param.CI = 0.95)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0298507 -0.0045290 0.0006747 0.0073960 0.0295073
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.953
## lower bound : 0.000, p = < 2.22e-16
## upper bound : 1.000, p = 0.00068248
## 95.0% CI : (0.885, 0.988)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.154229 0.059419 2.5956 0.01036 *
## log10(rootmass_av) 0.654704 0.042908 15.2584 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01015 on 153 degrees of freedom
## Multiple R-squared: 0.6034, Adjusted R-squared: 0.6008
## F-statistic: 232.8 on 1 and 153 DF, p-value: < 2.2e-16
# Quick plot of pgls
plot(log10(cor_av) ~ log10(rootmass_av), data = av.sal)
abline(pgls_cdmass)
# Plot average cornea diameter (mm) against cube root of mass (g)
plot_pgls_cdmass <- ggplot(av.sal, aes(y=cor_av, x=rootmass_av, text=Genus_Species)) +
geom_point(alpha = 0.9, aes(color = Adult_Habitat), size = 5, na.rm = TRUE) +
scale_color_manual(values = col_hab, name = "Adult habitat") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_text(size = 26),
axis.text.y = element_text(size = 26),
axis.title.x = element_text(size = 30, color = "black", margin = margin(t = 10)),
axis.title.y = element_text(size = 30, color = "black", margin = margin(r = 18)),
axis.ticks.length = unit(0.3, "cm"),
panel.border = element_rect(size = 2)
) +
scale_y_log10(name = "Cornea diameter (mm)") +
scale_x_log10(name = "Cube root of mass (g)", breaks = c(1, 3, 5), limits = c(0.5, 5.5)) +
geom_abline(slope = coef(pgls_cdmass)[[2]], intercept = coef(pgls_cdmass)[[1]], linetype = "solid", size = 1.5)
plot(plot_pgls_cdmass)
# Save the plot as a rectangle with specified width and height
ggsave("../Plots/plot_pgls_cdmass.png", plot = plot_pgls_cdmass, width = 9, height = 7.5, units = "in")
# PGLS cornea diameter vs. SVL
pgls_cdsvl <- pgls(log10(cor_av) ~ log10(svl_av),
data = sal.comp,
lambda = "ML", #uses Maximum Likelihood estimate of lambda
param.CI = 0.95)
# Diagnostic plots
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_cdsvl)
par(mfrow = c(1, 1))
# Likelihood plot for Pagel's lambda. Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval
lambda.cdsvl <- pgls.profile(pgls_cdsvl, "lambda")
plot(lambda.cdsvl)
# Print model output
summary(pgls_cdsvl)
##
## Call:
## pgls(formula = log10(cor_av) ~ log10(svl_av), data = sal.comp,
## lambda = "ML", param.CI = 0.95)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.034536 -0.004406 0.000953 0.007409 0.029583
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.958
## lower bound : 0.000, p = < 2.22e-16
## upper bound : 1.000, p = 8.8349e-06
## 95.0% CI : (0.904, 0.985)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.809986 0.104842 -7.7258 1.363e-12 ***
## log10(svl_av) 0.607606 0.042201 14.3980 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01066 on 153 degrees of freedom
## Multiple R-squared: 0.5754, Adjusted R-squared: 0.5726
## F-statistic: 207.3 on 1 and 153 DF, p-value: < 2.2e-16
# Quick plot of pgls
plot(log10(cor_av) ~ log10(svl_av), data = av.sal)
abline(pgls_cdsvl)
# Plot average cornea diameter (mm) against snout-vent length (mm)
plot_pgls_cdsvl <- ggplot(av.sal, aes(y=cor_av, x=svl_av, text=Genus_Species)) +
geom_point(alpha = 0.9, aes(color = Adult_Habitat), size = 5, na.rm = TRUE) +
scale_color_manual(values = col_hab, name = "Adult habitat") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_text(size = 26),
axis.text.y = element_text(size = 26),
axis.title.x = element_text(size = 30, color = "black", margin = margin(t = 10)),
axis.title.y = element_text(size = 30, color = "black", margin = margin(r = 18)),
axis.ticks.length = unit(0.3, "cm"),
panel.border = element_rect(size = 2)
) +
scale_y_log10(name = "Cornea diameter (mm)") +
scale_x_log10(name = "Snout-vent length (mm)", breaks = c(30, 70, 150), limits = c(25,160)) +
geom_abline(slope = coef(pgls_cdsvl)[[2]], intercept = coef(pgls_cdsvl)[[1]], linetype = "solid", size = 1.5)
plot(plot_pgls_cdsvl)
# Save the plot as a rectangle with specified width and height
ggsave("../Plots/plot_pgls_cdsvl.png", plot = plot_pgls_cdsvl, width = 9, height = 7.5, units = "in")
# PGLS svl vs. cube root of mass
pgls_svlmass <- pgls(log10(svl_av) ~ log10(rootmass_av),
data = sal.comp,
lambda = "ML", #uses Maximum Likelihood estimate of lambda
param.CI = 0.95)
# Diagnostic plots
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_svlmass)
par(mfrow = c(1, 1))
# Likelihood plot for Pagel's lambda. Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval
lambda.svlmass <- pgls.profile(pgls_svlmass, "lambda")
plot(lambda.svlmass)
# Print model output
summary(pgls_svlmass)
##
## Call:
## pgls(formula = log10(svl_av) ~ log10(rootmass_av), data = sal.comp,
## lambda = "ML", param.CI = 0.95)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0173001 -0.0037582 -0.0003353 0.0026855 0.0247753
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.921
## lower bound : 0.000, p = < 2.22e-16
## upper bound : 1.000, p = 1.8264e-12
## 95.0% CI : (0.847, 0.963)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.620475 0.033299 48.664 < 2.2e-16 ***
## log10(rootmass_av) 0.996038 0.026331 37.827 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.005711 on 153 degrees of freedom
## Multiple R-squared: 0.9034, Adjusted R-squared: 0.9028
## F-statistic: 1431 on 1 and 153 DF, p-value: < 2.2e-16
# Quick plot of pgls
plot(log10(svl_av) ~ log10(rootmass_av), data = av.sal)
abline(pgls_svlmass)
# Plot average snout-vent length (mm) against cube root of mass (g)
plot_pgls_svlmass <- ggplot(av.sal, aes(y=svl_av, x=rootmass_av, text=Genus_Species)) +
geom_point(alpha = 0.9, aes(color = Adult_Habitat), size = 5, na.rm = TRUE) +
scale_color_manual(values = col_hab, name = "Adult habitat") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_text(size = 26.3),
axis.text.y = element_text(size = 26.3),
axis.title.x = element_text(size = 30.3, color = "black", margin = margin(t = 10)),
axis.title.y = element_text(size = 30.3, color = "black", margin = margin(r = 12)),
axis.ticks.length = unit(0.3, "cm"),
panel.border = element_rect(size = 2)
) +
scale_y_log10(name = "Snout-vent length (mm)", breaks = c(30, 70, 150), limits = c(25,160)) +
scale_x_log10(name = "Cube root of mass (g)", breaks = c(1, 3, 5), limits = c(0.5,5.5)) +
geom_abline(slope = coef(pgls_svlmass)[[2]], intercept = coef(pgls_svlmass)[[1]], linetype = "solid", size = 1.5)
plot(plot_pgls_svlmass)
# Save the plot as a rectangle with specified width and height
ggsave("../Plots/plot_pgls_svlmass.png", plot = plot_pgls_svlmass, width = 9.4, height = 7.5, units = "in")